% Approximate Hellinger distance in bivariate example under shock-sign
% restriction.

clear variables
close all
clc

curFolder = pwd;

rng(23032021); % Set random seed

T = 3; % Length of time series
n = 1000; % Number of gridpoints over which to evaluate likelihood
k = 1; % Period in which first structural shock is positive
K = 1e6; % Monte Carlo sample size for approximating Hellinger distance

% Set true values of structural parameters.
A0 = [1 0.5; 0.2 1.2]; % Implies sig21 < 0

% Compute reduced-form parameters
% Lower-triangular Cholesky factorisation of innovation covariance matrix.
Sigmatr = chol(inv(A0)*inv(A0)','lower'); 
Sigma = Sigmatr*Sigmatr'; % Covariance matrix of VAR innovations.
sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

% True value of theta in orthogonal reduced form.
thetaTrue = asin(A0(1,2)*sig22); % sig21 < 0

% Construct grid for theta.
thetaGrid = linspace(-pi+eps,pi-eps,n)';

%% Approximate Hellinger distance for conditional likelihood.

hellingerC = zeros(n,1);

for jj = 1:K % For each Monte Carlo sample

    % Draw a time series of structural shocks such that the first 
    % structural shock in the kth period is positive.
    epst = randn([T,2]);
    flag = 0;

    while flag == 0

        epst(k,1) = randn;
        flag = (epst(k,1) >= 0);

    end

    % Use shocks and structural parameters to generate data.
    yt = (A0\epst')';

    h = sig21*yt(k,1)-sig11*yt(k,2);
    C = sig22*yt(k,1)/h;

    if h < 0 % Case (I)

        theta_cis = [atan(max(sig22/sig21,C)), pi + atan(min(sig22/sig21,C))];

    elseif h > 0 && sig22/sig21 < C % Case (II)

        theta_cis = [atan(sig22/sig21), atan(C)];

    elseif h > 0 && sig22/sig21 > C % Case (III)

        theta_cis = [pi + atan(C), pi + atan(sig22/sig21)];

    end

    % Add contribution of draw to Hellinger distance.
    hellingerC = hellingerC + ...
        ((theta_cis(1) <= thetaGrid) & (thetaGrid <= theta_cis(2)))./K;

end
  
HDCsq = 2*(1 - hellingerC); % Squared Hellinger distance

%% Approximate Hellinger distance for unconditional likelihood.

hellinger = zeros(n,1);

for kk = 1:K
    
    % Draw a time series of structural shocks and check whether shock-sign
    % restriction holds.
    epst = randn([T,2]);
    flagSR = (epst(k,1) >= 0);

    % Use shocks and structural parameters to generate data.
    yt = (A0\epst')';
   
    unLik = zeros(n,1);
    
    for ii = 1:n % For each value of theta

        theta = thetaGrid(ii);

        % Value of first structural shock given value of theta and 
        % realisation of data in period k.
        eps1k = (1/(sig11*sig22))*(sig22*yt(k,1)*cos(theta) ...
            + (sig11*yt(k,2)-sig21*yt(k,1))*sin(theta));

        % Check whether narrative restriction and sign normalisations are 
        % satisfied at value of theta.

        if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && (eps1k*flagSR -eps1k*(1-flagSR) >= 0)
            
                % Add contribution to Hellinger distance.
                hellinger(ii) = hellinger(ii) + 1/K;

        elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
                && (eps1k*flagSR -eps1k*(1-flagSR) >= 0)

                % Add contribution to Hellinger distance.
                hellinger(ii) = hellinger(ii) + 1/K;

        end

    end   
    
end

HDsq = 2*(1 - hellinger);

figure;
h1 = plot(thetaGrid,HDsq,'color',[255 0 0]./255,'LineWidth',2);
hold on;
h2 = plot(thetaGrid,HDCsq,'color',[255 0 255]./255,'LineWidth',2,...
    'LineStyle','--');
h3 = line([thetaTrue thetaTrue],[0 2],'color',[76 0 153]./255,'LineWidth',2);
xlim([-pi pi]);
xticks([-pi -pi/2 0 pi/2 pi]);
xticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'});
xlabel('\theta');
ylim([0 2]);
title('Shock-sign Restriction');
Ax = gca;
Ax.FontSize = 14;
cd('Figures');
print('ShockSignHellinger','-depsc');
print('ShockSignHellinger','-dpng');
cd(curFolder);  

cd('Results');
save('ShockSignHellinger_results.mat');
cd ..